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First-order phase transition and the equation of state in a 2D granular fluid 
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We present experimental evidence for a first-order freezing/melting phase transition in a nonequi- 
librium system — an oscillated two-dimensional isobaric granular fluid. The steady-state transition 
occurs between a gas and a crystal and is characterized by a discontinuous change in both den- 
sity and temperature. It is suppressed if the number of particles is incommensurate with the cell 
size, shows rate-dependent hysteresis, and obeys the Lindemann criterion for melting. Further, the 
measured equation of state both above and below the phase transition compares well with theory. 

PACS numbers: 45.70.-n, 51.30.+J, 51.10.+y, 64.70.Hz 



Introduction: Granular materials are fascinating sys- 
tems, which have tremendous technological importance 
and numerous applications to natural systems. They also 
represent a serious challenge to statistical physics as an 
extreme example of a system far from equilibrium. As 
for any macroscopic system, the total number of modes 
in a collection of grains is on the order of Avogadro's 
number Na- However, in granular systems a very small 
number of modes, specifically the translational and ro- 
tational mode of the N macroscopic grains which make 
up the system, can be preferentially excited by exter- 
nal forces like moving walls. The energy in these ~ N 
modes can be many orders of magnitude above the av- 
erage energy per mode (i.e., the temperature) for the 
remaining Na modes. The study of the relaxation of 
this extremely far-from-equilibrium system toward equi- 
librium is the challenge of granular statistical mechanics. 
In this paper, we experimentally examine the nonequilib- 
rium steady-state (NESS) created by a balance between 
relaxation (dissipation) and injection of energy (heating), 
and specifically, a first-order melting/freezing phase tran- 
sition and the equation of state in a two-dimensional (2D) 
inelastic hard sphere system heated (excited) from below 
under isobaric conditions. 

The elastic hard sphere system is the simplest system 
which undergoes a first-order phase transitionpJQ- This 
transition has been seen in simulation Q and experi- 
ments in colloids 4J. However, in these systems energy is 
conserved, and the concept of free energy is well defined. 
In dissipative systems, such as inelastic hard spheres or 
granular systems, the concept of free energy is not es- 
tablished |5j, even in a NESS, since there is a constant 
flow of energy through the system. The idea of applying 
thermodynamics to a NESS is at the forefront of cur- 
rent statistical physics Q . The demonstration of exper- 
imental systems which undergo first-order phase transi- 
tions under NESS conditions points to the universality of 
entropy-production and free-energy concepts even in the 
absence of energy conservation. 

Experiment: We place N (26-85) spherical stainless 
steel ball bearings of diameter D = 3.175 mm in a con- 
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FIG. 1: (a)-(c) Photographs of a 2D granular layer, with freely 
floating (isobaric) weight: (a) R = 3.5, T = 3.88, Tvv = 0.12, 
disordered with crystalline regions, (b) R = 4, V = 7.90, 
Tvv = 0.12, crystal, (c) R = 4, T = 8.10, Tvv = 0.59, gas. 
(d) The radial distribution function for states (a)-(c). The 
line style is indicated below each image. The vertical lines 
show next-nearest-neighbor distances. The short line is \/3~D 
and the long line is \/3 times the average nearest-neighbor 
distance in (b). The inset shows the first peak at full scale. 



tainer 17.5 D wide by 20 D tall by 1 D deep as shown in 
Fig. QJa-c). We define the number of rows R = N/17, 
where 17 is the number of particles to fill an entire row. 
A thin plunger slides through a slot in the bottom of the 
cell and oscillates sinusoidally to excite (heat) the par- 
ticles from below. The driving is characterized by the 
nondimensional maximum acceleration Y — A(2irf) 2 /g, 
where A is the maximum amplitude of the plunger, / 
is the frequency and g is the acceleration of gravity. A 
freely floating weight confines the particles from the top, 
allowing the volume to fluctuate but providing constant 
pressure conditions. When Y is large the granular system 
acts like a gas, as shown in Fig. [He). When Y is small 
and R is an integer a crystalline state develops [see Fig. 



IHb)]. Using high-speed digital photography we measure 
the positions of the plunger, the weight, and all of the 
particles in the cell with a relative accuracy of 0.2% of 
D or approximately 6/xm at a rate of 840 Hz. We track 
the particles from frame to frame and assign a velocity 
to each one, typically ~ D/5 per frame. 

From these data we measure the volume, density, and 
granular temperature or average kinetic energy per par- 
ticle. In a normal gas, the temperature is isotropic 
and the kinetic energy in vertical velocity is the same 
as that in horizontal velocity. In a granular system, 
due to dissipation, at least two temperatures are needed 
in two dimensions, a horizontal temperature Thh = 
l/2m((v x — (v x )) 2 } and a vertical temperature Tyv — 
l/2m((v y — (v v )) 2 ), where m is the mass of the particles 
and v x and v y are the horizontal and vertical velocities 
respectively. We want to assign a single value to each 
temperature, so we take the averages over all particles 
and time (~ 50 periods). The most obvious implemen- 
tation of the averages (...) is in the inertial lab frame; 
however, this produces a finite temperature even if the 
plunger, particles, and weight all move with the same ve- 
locity. Taking the averages in the center of mass frame of 
the particles solves this problem, but the velocities are no 
longer simply related to the collision velocities with the 
weight used in the derivation of the equation of state. We 
have chosen to evaluate the averages in the frame of the 
weight. This choice solves both problems. The ultimate 
solution is to allow the temperature, density, and pres- 
sure to vary in space and time, and look at the space and 
time dependent equation of state. However, our current 
approach is much simpler, and the experimental equation 
of state compares well with standard kinetic theory 0. 

We measure the average density, horizontal and verti- 
cal temperatures, and radial distribution functions g(r) 
for three different constant pressures, for R from 1.5 - 
5, and for Y from 0-30. To prepare the system initially 
in the densest state, we increase Y to 30 and then slowly 
lower the acceleration to zero. Then we alternate taking 
1024 pictures with a small increase in T and a variable 
time delay (typically one second) until T = 30. Then 
we repeat the same process while decreasing T to zero. 
Typical data runs are shown in Fig. 

Results: The behavior of the system is quite different 
depending on the number of particles in the cell. If R is 
an integer then we observe a first-order hysteretic phase 
transition as we change T. Figure E^a) shows a typical 
example for R = 4. Below r = 1, nothing happens - 
the temperatures are zero and the 2D volume fraction v 
is constant. As we increase Y for 1 < Y < 8, Tyv rises 
linearly. Initially, Thh — Tvv/% but around r = 3 the 
ratio drops and stays around 20%, as shown in Fig. Hfc). 
v remains nearly constant for T < 2 and then falls linearly 
until r = 8. At r = 8, there is an abrupt change in all of 
the measured quantities. Tyv and Thh rise by factors of 
5 and 10 respectively, and the ratio reaches nearly 50%. 
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FIG. 2: (a)-(b) Plots of the volume fraction (thick) and aver- 
age vertical (dashed) and horizontal (thin) temperatures as a 
function of increasing and decreasing P at 50 Hz, under iso- 
baric conditions, (a) First-order phase transition for R = 4, 
showing a discontinuities in density and temperature as well 
as hysteresis, (b) For R — 3.5, the curves are smooth and non- 
hysteretic. (c) Plot of the horizontal to vertical temperature 
ratio for R — 3.5 (dashed) and R = 4 (solid). 



v drops by more than 15%. For T > 8, Tyy and Thh 
continue to rise, the ratio reaches 70% by Y — 30 (not 
shown), and v continues to drop. As Y is lowered, the 
process is reversed, but the transition point shows 25% 
hysteresis, occurring at Y = 6.1. The hysteresis is rate- 
dependent, and if the rate of change in Y is increased by 
a factor of 100 the hysteresis is lost. Similar behavior is 
seen for R = 2, 3, and 5, although the size of the jumps 
decrease as R is decreased. For R = 3.5 the situation is 
quite different, as shown in Fig. |2Jb). All of the measured 
quantities change continuously as Y is changed, and no 
hysteresis is observed. The ratio of Tyv to Thh stays 
near 50% until Y = 1. Similar behavior is seen for R = 
1.5, 2.5, and 4.5. 

g(r) is shown in Fig. ^d) for the R = 4 states just be- 
fore (r = 7.90) and just after (r = 8.10) the transition. 
While r only changes by 2.5% g(r) changes dramatically. 
g(r) for the gas state is typical for a hard sphere fluid near 
a freezing transition, including the broad peak just below 
r = 2D that is a precursor for freezing [8j. g(r) for the 
crystalline state is typical of an expanded hard sphere 
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crystal at nonzero temperature, showing sharp peaks at 
r = D, 2D, and 3D. There is a big peak at the next- 
nearest neighbor distance r = y/3D with a splitting due 
to the fact that the average nearest neighbor distance 
is greater than D at a nonzero temperature. g(r) for 
R = 3.5 at r = 8 is not shown but is clearly fluid-like. 
Even at r = 3.88 where the Tyy is equal to the crys- 
tal melting temperature, g(r) is intermediate between a 
crystal and fluid, with small peaks at integer separations, 
as shown in Fig. dd). In the image [Fig. Wi a )] crystalline 
regions can be seen, but there are gaps, holes, and dis- 
locations. Information from g(r) is summarized in Fig. 
EJb-c). In all of the crystallizing systems particles are ex- 
cluded from a band at r = 1.5D ± 0.1D and the average 
number of next nearest neighbors is near the maximum 
as shown in Fig. EJb-c) . The integer system is clearly 
in a crystalline state, but it is not clear if the half-integer 
system ever reaches a crystalline state. Another clas- 
sic measure of melting is the Lindemann criterion || , in 
which the average root-mean-squared particle displace- 
ments 7„j = y{\f— WFF)/D > lc- lm is plotted against 
T in Fig. C^a). We find a value of j c between 0.10-0.12 
is consistent with all of the integer data. This also pre- 
dicts that the half-integer systems freeze, but at a lower 
r, around 4 for R — 3.5. However, the predicted freez- 
ing temperatures for all data fall in a band between the 
crystal freezing and melting temperatures (Tyy — 0.1 - 
0.6). While the Lindemann criterion suggests that all of 
the systems freeze, it says nothing definitive about the 
structure of the solid and the radial distribution func- 
tions suggest that freezing is a continuous process for the 
half-integer systems. Recent work on elastic hard spheres 
using density-functional theory by Both and Hong [Toj j 
predicts a critical melting temperature determined by a 
constant /^o- From our melting and freezing temperature 
we find a range of /iq between 10-60 which is consistent 
with their findings. 

To further exploit the analogy to thermodynamics, 
we have examined the average equation of state for 
this system. From Tyy, the average number density 
n = Ah>/(jtD 2 ), and average pressure determined from 
the mass of the weight and half of the mass of the parti- 
cles, we experimentally measure the compressibility fac- 
tor y from the equation of state for this 2D granular 
fluidtH, P = nT vv (l + x (y))\ X = P/(nT vv ) - 1. X 
can be theoretically related to g(D), through the viral 
for the pressure [l2j. In an inelastic hard sphere fluid 
X{v) = otvg(D;v) = aGiv), where a = (1 + e) and e is 
the coefficient of restitution. For an elastic hard sphere 
fluid e = 1 and a = 2. A number of forms for G are 
available in the literature (see [13J for a recent discus- 
sion). All are similar and we use the simple Gt developed 
by Torquato |14| . which is an analytical fit to molecular 
dynamics simulation at high v and the Carnahan and 
Starling |15| geometric series approximation to the first 
few viral coefficients at low v. 
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FIG. 3: Plots of the (a) Lindemann ratio and the average 
number of particles separated by (b) %/3D ± 0.1D and (c) 
1.5D±0.1L> versus r for R = 4 (thick) and 7? = 3.5 (dashed). 



A plot of x versus volume fraction is shown in Fig. 0]for 
R = 3.5 and 4. x remains ~ 10% of 2Gt up to v = 0.55, 
R = 3.5 and 4. Above 0.55 the half-integer system drops 
~ 10% below 2Gt, but both X an d its derivative are 
continuous functions of v up to x — 200. Above this 
value, r is near 1, the temperatures are almost zero, and 
X strongly diverges. This behavior is expected, since our 
equation of state completely neglects the elastic forces 
between the particles. When the temperature goes to 
zero the weight is held up by elastic energy, not kinetic 
energy. However, the point at which this occurs is a very 
high density in which the pressure due to the density of 
the gas through the compressibility factor is 200 times 
the value of a dilute gas. The behavior above v = 0.6 
is quite different for the integer systems. As v increases 
there is a jump in both x ancl v signaling the phase tran- 
sition. Above this volume fraction the system is in an 
expanded crystalline state like that shown in Fig. db). 
The size and position of the gap depend on a number of 
parameters, including the history, number of rows, and 
confining pressure. In the crystalline state x increases 
until r is nearly one and the maximum recorded packing 
fraction is reached. At this point x diverges strongly as 
elastic forces become important. As with the half-integer 
case, the granular temperature must be very low and the 
density must be very close to the maximum before elastic 
consideration are important. This suggests that contin- 
uum theories based solely on kinetic considerations may 
have a range of applicability up to very high density and 
very low temperatures. 
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FIG. 4: Log-linear plot of the compressibility factor as a func- 
tion of volume fraction, for R = 3.5 (dashed), R — 4 (thick), 
and 2Gt (thin, see text). The dashed (R = 3.5) and thick 
(7? = 4) curves in the inset show an experimental measure- 
ment of a e xp (see text) as a function of volume fraction, a — 2 
for elastic hard spheres. The large dots (R = 4) and small 
dots (R = 3.5) show ctr using Gt- 



To see the comparison between \ and Gt more clearly, 
the inset of Fig. 0] plots x/Gt which should be 2 for elas- 
tic hard spheres and less for inelastic hard spheres. We 
can also measure G from the experimental g(D;v) and 
compare that to the experimental value of \. G exp deter- 
mined in this manner is only accurate for v < 0.8; above 
this value the peak near D begins to split due to small 
imperfections in the crystal lattice, and the peak height 
stops growing. From G exp we obtain a fully experimental 
measure of a. The result is also shown on the inset of 
Fig. 0] a exp is consistent with our direct measures of e 
from particle tracking. We find values of e over a broad 
ranging from 0.3 to 1.3, and a general trend of lower val- 
ues for higher collision velocities. Values of e > 1 are 
possible since we cannot measure particle spin. 

Discussion: The first-order phase transition seen in the 
integer number of rows is qualitatively different from that 
of an elastic hard sphere system. In such a system, just 
as in an ordinary gas, at the phase transition there is a 
discontinuous change in the density, but the temperature 
would be unchanged. This is a unique feature in granular 
systems since it suggests that there is the possibility of 
steady states in which two phases co-exist, but at differ- 
ent nonzero temperatures. States like these have been 
seen experimentally and in simulation |l7| . This co- 



existence does not violate thermodynamics since a steady 
state is not the same as equilibrium, but it can have pro- 
found effects on continuum theories since it predicts dis- 
continuities in the temperature field could exist in steady 
state without the need for energy or particle flow. That 
is, the system is mechanically stable since there are two 
states with the same pressure but different densities and 
internal energies. This would also interfere with the stan- 
dard continuum assumption of near Maxwell-Boltzmann 
distributions, since near the gas/crystal border two dif- 
ferent velocity distributions can co-exist. Any spatial av- 
erage would include some of each distribution. So while 
the equation of state data suggests that theories based 
solely on kinetic consideration are applicable on either 
side of the phase transition, some modification will be 
needed to account for the discontinuities created by phase 
transitions. 
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